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Abstract 

A new molecular dynamics model in which the point charges on atomic sites are 
allowed to fluctuate in response to the environment is developed and applied to water. 
The idea for treating charges as variables is based on the concept of electronegativity 
equalization according to which: (a) The electronegativity of an atomic site is depen- 
dent on the atom's type and charge and is perturbed by the electrostatic potential it 
experiences from its neighbors and (b) Charge is transferred between atomic sites in 
such a way that electronegativities are equalized. The charges are treated as dynamical 
variables using an extended Lagrangian method in which the charges are given a ficti- 
tious mass, velocities and kinetic energy and then propagated according to Newtonian 
mechanics along with the atomic degrees of freedom. Models for water with fluctu- 
ating charges are developed using the geometries of two common fixed-charge water 
potentials: the simple point charge (SPC) and the 4-point transferable intermolecular 
potential (TIP4P). Both fluctuating charge models give accurate predictions for gas- 
phase and liquid state properties, including radial distribution functions, the dielectric 
constant, and the diffusion constant. The method does not introduce any new inter- 
molecular interactions beyond those already present in the flxed charge models and 
increases the computer time by only a factor of 1.1, making this method tractable for 
large systems. 
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1 Introduction 



In simple molecular force fields the intramolecular electronic structure is often modeled by 
point charges fixed on well defined sites in the molecular frame. The charges are constant and 
thus cannot change in response to changing electrostatic fields which arise from movement 
of the atoms during the simulation. In reality, molecular electronic structure can be strongly 
influenced by the molecular environment. For example the total dipole moment of water 
changes from 1.85 Debye in the gas phase to approximately 2.5 Debye in the liquid phase. 
Thus the charges used in simulations based on fixed charge force fields must reflect average 
or mean field charge values for the particular phase and are not transferable to different 
thermodynamic states or to different media. In addition, the self-energy involved in the 
change in charge accompanying the transition from gas phase to liquid phase is commonly 
neglected in most force field parametrizations. This "missing term" in fixed charge pair 
potentials can be significant (2 to 5 kcal/mole for water) [?]. Charge induction effects are 
not pair-wise additive and improved models must go beyond pair potentials. The purpose 
of this paper is to present a new simulation method in which the charges are responsive to 
environmental changes. 

The approach taken here combines the electronegativity equalization (EE) method for 
solving for atomic charges and the extended Lagrangian method for treating fictitious degrees 
of freedom as dynamical variables. The calculation of atomic electronegativities using density 
functional theory, basis set methods, or empirical data and the use of this information to 
estimate charges for large molecules has a long history[?, ?, ?, ?, ?, ?]. The electronegativity 
of an atomic site is dependent on its charge and the electronegativities of the neighboring 
atoms. Parr has shown that the MuUiken electronegativity (xi) of an isolated atom i is the 
negative of the chemical potential (/ij) of the electron gas surrounding its nucleus. 



where E is the ground state energy, N is the number of electrons in the atom (treated as 
a continuous variable), Q is the charge on the atom, and e is the elementary charge. Use 



dE _ 
dN ~ 



dQi 



dE 



(1.1) 
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has been made of the fact that Q is related to iV by Q = —e{N — Z) where Z is the 
atomic number of the atom. In a many atom system, the electron gas will equilibrate with 
the instantaneous positions of the nuclei in such a way that the electrochemical potential 
of the electron gas will be equal at all atomic sites. In this picture, electrons will then 
move among atoms from regions of low electronegativity (or high electrochemical potential) 
to regions of high electronegativity (low electrochemical potential). For the ground state 
electronic configuration, the electrochemical potentials are equal. The approach we take is to 
treat charges on the molecular sites as dynamical variables by introducing fictitious kinetic 
energy terms and self energy terms for these charges into the Lagrangian for the system 
along with Lagrange constraints representing various conditions of electroneutrality. In this 
extended Lagrangian approach [?, ?, ?, ?] the charges are propagated according to Newtonian 
mechanics in a similar way to the atomic degrees of freedom. Although the fluctuating 
charge model provides a convenient starting point for the discussion of complex solutes like 
n-methylacetamide and proteins, we focus here on its application to the simulation of neat 
water. 

Liquid water was chosen for the first application of this fluctuating charge (or fluc-q) 
model because charge polarization effects should be important for water. Simple water 
models, such as the simple point charge (SPC)[?] or the 4-point transferable intermolecu- 
lar potential TIP4P[?], with a Lennard- Jones interaction between oxygen atoms and three 
charge sites with flxed liquid state charges, can give accurate predictions for many equi- 
librium properties of liquid water, including the energy and radial distribution functions. 
Due to their simplicity and relative accuracy, these are perhaps the two most widely used 
water potentials. However, the translational and rotational time scales of these models are 
too fast, although the SPC/E[?] reparameterization gives improved relaxation times. More 
importantly, these models are unable to deal accurately with heterogeneous environments. 
Simulations of sodium octanoate micelles in SPC water predict too much penetration of wa- 
ter molecules into the micelle [?, ?] in contrast to experiments [?]. The purpose of this paper 
is to add fluctuating charges to the SPC and TIP4P potentials and to thus devise a model 
which has improved static and dynamical properties relative to the fixed charge models and 
which should is easily extendible to more complex solutions. 



3 



Electrical induction can also be described to lowest order using fixed gas phase charges 
and point polarizabilities. Many dipole polarization models have been used to simulate 
liquid water [?, ?, ?, ?, ?, ?, ?, ?, ?]. The fiuc-q model presented here is an alternative 
which differs from the dipole polarizable models in two respects. First, the fluc-q models 
have polarizabilities to all orders in the charge moments and not only dipolar polarizability. 
In addition, the dipole polarizability models introduce a new interaction (the 1/r'^ dipole- 
dipole interaction) and in order to solve for the induced dipole moments, either the induced 
dipole equations are solved iteratively, by matrix inversion, or the polarizations are treated 
in an extended Lagrangian framework. The iterative solution method increases the cost by 
a factor of two[?] and the extended Lagrangian methods by a factor of two[?] to four[?]. The 
fluc-q method introduces no new intermolecular interactions beyond the fixed-charge models 
and increases the cpu time by a factor of only 1.1. Of course the dipole polarizability model 
can also be cast in terms of Drude dispersion oscillators, a system that also lends itself to 
treatment by the extended Lagrangian method. 

This paper is organized as follows as follows. Section ^ describes the fluc-q method and the 
form of the water-water interactions. Section ^ describes our implementation of molecular 
dynamics. Section ^ describes the results of these models and Section |^ summarizes the 
conclusions. 

2 Dynamical Fluctuating Charge Models 

The central idea for treating charges as dynamical variables is based on the electronegativities 
of atomic sites. Parr has shown, using Kohn-Sham theory, that in an atom, the atomic 
electrons, regarded as an electron gas, have a chemical potential which is the negative of the 
MuUiken electronegativity [?]. In a many-atom system the full electron gas will distribute 
itself so that its electrochemical potential takes the same value at every nuclear site. This 
principle of electronegativity equalization (EE) was first proposed by Sanderson[?]. If a given 
site moves so that it feels a different electrostatic potential it will take on a different charge. 
In this way the charges on molecular sites will respond to the environment. 
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In the isolated atom the energy of creating a partial charge, Qa, can be expanded to 
second order as 

E(Q^) = E^(0) + x^Q^ + ij^^Qi (2.1) 

where and J^^, are parameters dependent on the atom type. Values of Xa J^a 
be calculated using basis set, density functional theory methods or empirical data: Xa 
the MuUiken electronegativity (per electronic charge e) and J^^ is twice the hardness of the 
electronegativity of the isolated atom. The energy of a system of Nmoiec molecules each with 
Natom atoms is 

^molec Natom / 1 \ 

U({Q},{r})= E (Ea(0)+£Q,a + -jLQL)+E Ja/3(w)QiaQ.73+E V(W) 

1=1 a=l ^ ^ ia<jl3 ia<jl3 

(2.2) 

where Eq,(0) is the ground state energy of atom a, Tiajp is the distance, Ja/sijiajp) the 
Coulomb interaction and V(rjQ,j^) is any additional non-Coulombic interaction between ia 
and jp. The electronegativity per unit charge of atom A is given by 

The charges, by the EE principle, are then those for which the electronegativities are equal. 
This is equivalent to minimizing the energy, subject to a charge neutrality constraint. Since 
the potential is quadratic in the charges, the minimization will lead to a set of coupled linear 
equations for the charge. 

The charges are not independent variables since there is a charge conservation constraint. 
For uncharged molecular systems, the constraint can be of two types: 

1. The entire system is constrained to be neutral, so individual molecules can carry a 
non-zero charge because there can be intermolecular charge transfer, 

^molec Natom 

E E = 0- (2-4-a) 

1=1 a=l 
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2. Each molecule is constrained to be neutral, so there is no intermolecular charge transfer, 

^ atom 

E Q- = 0- (2-4-b) 

a=l 

With intermolecular charge transfer, the chemical potentials of all the atoms of the system 
will be equal. Without charge transfer, the chemical potential of an atom will only be equal 
to the chemical potential of atoms on the same molecule. The simplest way to treat the 
charge neutrality constraint is to treat the charges as independent and use the method of 
undetermined multipliers to enforce the constraint. The Lagrangians for cases 1 and 2 are: 

Li= ^rnJl+ E E 2MQQ'.-U({Q},{r})-A E E (2-5-a) 

i=l a=l i=l a=l i=l a=l 



and 



^2= E E 2^arL+ E E 2MQQL-U({Q},{r})- E E Q- (2-5-b) 

i=l a=l 1=1 a=l 1=1 a=l 

where m^ is the mass of atom a and Mq is a fictitious charge "mass", which has units of 
energy-sec^ /charge^ and the A's are Lagrange multipliers. The nuclear degrees of freedom 
evolve according to Newton's equation 



and the set of charges evolve in time according to 

. . ^ c}U({Q},{r}) . ~ . ..7, 

MqQia = Ai = -Xia - Ai (2.7) 

where Aj is the Lagrange multiplier for the charge neutrality constraint, given either by 
Eq. ( |2.4-a| ) or Eq. ( p.4-b| ). It should be noted that if the total charge in the liquid is a 



constant of the motion, then 



E E Q- = (2-8-a) 

1=1 a=l 
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whereas if the total charge on molecule i is a constant of the motion, then for all i 

^ atom 

E Q- = 0. (2.8-b) 
Substitution of Eq. ( |2.7| ) into each of the above two equations yields, respectively, 

^ = ^ E E X^. (2.9-a) 

-'■ ^ molec-'- ^ atom {—i a=l 

where A is equal to the negative of the average of the system's total electronegativity, and 

Natom 

K = E X^a (2.9-b) 

-^^ atom a=l 

where Aj is the negative of the average electronegativity on molecule i. Substitution of 
Eqs. ( |2.9-a| ) and ( |2.9-b| ) into the equations of the motion for the charges gives: 



MQQ^a = -- E E iX^a-X^p) (2.10-a) 

i=l (3=1 



-^molec-^atom 



and 



MqQ,^ = -- E iX^a-X^p)■ (2.10-b) 

13=1 



^atom 



Whether or not charges are allowed to transfer between molecules or just between atoms 
on the same molecule makes little difference in the algorithm. The force on the charge is 
simply the difference between the average electronegativity and the electronegativity at that 
site. For example, if the electronegativity is greater than the average, then the force acts 
to decease the charge until the electronegativities are all equal. In the present application, 
we have included a charge neutrality constraint on each water molecule, there is no charge 
transfer between molecules (case 2). 

There have been many different applications of the EE principle, with differing values 
for the parameters xa and Jaa (some with higher order terms) which have been applied to 
a variety of molecules[?, ?, ?, ?, ?]. Some of these applications have been used as potential 
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input into simulations. However, in order to fully treat the charge fluctuations which arise 
in the course of a simulation, each new conflguration would require the calculation of a new 
set of charges. This can be done simply using the extended Lagrangian method to derive the 
equations of motion for Hamiltonians which depend on auxiliary degrees of freedom [?, ?]. 

The above approach makes no allowance for the possibility that there can be barriers 
to the charge transfer preventing charge equilibration. Thus for example two atomic sites 
separated by a great distance will in reality not transfer charges because the probability of 
tunneling through a wide barrier is small, yet the above model will allow charge transfer. 
This poses a major problem in principle when molecules are separated in vacuo, a problem 
that also exists in the Car-Parrinello method[?]. For this reason, we have restricted the 
charge equilibration to intramolecular charge transfer. It would be useful to generalize this 
model to include such kinetic restrictions on charge transfer. 

The charge mass, Mq, a flctitious quantity, should be chosen to be small enough to 
guarantee that the charges readjust very rapidly to changes in the nuclear degrees of freedom. 
This is equivalent to the Born-Oppenheimer adiabatic separation between the electronic and 
nuclear degrees of freedom. When Mq is sufficiently small there will be essentially no thermal 
coupling between the nuclear and electronic degrees of freedom. For numerical convenience 
we choose the mass small enough to satisfy the foregoing requirement yet large enough so 
that the time-step required for the solution of the equations of motion is not too small. The 
charge degrees of freedom are to remain near zero Kelvin, since they are to be near the 
values which minimize the electrostatic energy. Within the context of molecular dynamics, 
this can be achieved by Nose thermostating[?] the charge at a much lower temperature 
(~5 K) than the nuclei as is done in Car-Parrinello ab initio molecular dynamics[?, ?, ?]. 
However, with a charge mass of 10.0 for the TIP4P-FQ model and 11.6 (psec/e)^ kcal/mole 
for the SPC-FQ model and a 1 femtosecond time step for both models, there is almost no 
thermal coupling and the charge degrees-of-freedom remain at a temperature of less than 6 
Kelvin for the duration of a 50 picosecond simulation. For simulations of this duration or 
shorter, no thermostating is needed to keep the charge degrees of freedom near zero Kelvin 
and the atomic degrees of freedom near the desired, much higher temperature. With this 
time step and Mq, the fluc-q models have the same energy conservation as the fixed charge 
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models. However, that small charge mass requires the use of a time step no greater than 
about a femtosecond. This problem could be surmounted using multiple time scale molecular 
dynamics [?]. 

Following Rappe and Goddard[?], the Coulomb interaction, Jjj(r), for intramolecular 
pairs is taken to be the Coulomb overlap integral between Slater orbitals centered on each 
atomic site, 

J,,(r) = / dvdv,\<PnX^^)\ \ _^ _ M n.ijl)?- (2.11) 

The Slater orbitals are given by 

0„^(r) = A,r"»-ie-^»'- (2.12) 

and are characterized by a principal quantum number, nj, and an exponent Q. Aj is a 
normalization factor. The value of Jji(r) for r=0 is J° and therefore the value of Q uniquely 
determines J°. For hydrogen, nj7=l and ^%h=\Ch and for oxygen, no=2 and Jqo^^Co- 
The intermolecular Coulomb interaction is set equal to the pure Coulomb interaction, 1/r, 
for consistency with other force fields. The interaction Jhh is shown in Figure |l|. 

Two different water geometries were used, corresponding to the commonly used SPC[?] 
and TIP4P[?] water models. Both of these models have 3 charged sites, two positive charged 
hydrogen sites and a negative charged site, and a Lennard- Jones interaction between oxygen 
sites, 



ULj(r) = 4e 



(T\i2 fa 



(2.13) 



. r / \ r , 

with a well depth, e, and a diameter, a. The SPC potential places the negatively charged site 
on the oxygen position; the TIP4P potential places this site (called the "M-site") a distance 
of 0.15 A from the oxygen position along the dipole direction toward the center-of-mass. The 
0-H bond length and H-O-H bond angle for the potentials are listed on Table |l]. The TIP4P 
model has the added complexity (and computational cost) of an additional interaction site, 
but has the correct water geometry. The potential energy contains the Lennard- Jones part 
(Eq. ( p.l3| )) and the electrostatic part (Eq. (p.2|)) and since we are defining the energies 
relative to the isolated gas-phase energy, the gas phase energy, E^p needs to be subtracted. 
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For the isolated gas-phase water molecule, the charge constraint gives Qo = —^Qh and it is 
straightforward to find that the charge which minimizes the energy is 



and the gas phase energy is thus 



(2.14) 



9P 



(2.15) 



where iur is the distance between the hydrogen and the M-site for the TIP4P model and 
for the SPC model, it is the distance between the hydrogen and the oxygen site. We have 
dropped the charge independent term (Ei(0)), taking this as our definition of the zero of 
energy. The total energy for ^moiec molecules is a sum of the Lennard- Jones part, the inter- 
molecular Coulomb part, an intramolecular self-energy and the gas- phase energy correction, 
to give 



1=1 j<ci \ 

N 7 



a 



JiO,jO 



12 



1 



JiOJO, 



3 3 

+ E E QiaQjis/'CiaJP 
a=l/3=l 



+ E E xlQia + 7^ E Q-Q 



^molec^gp 



(2.16) 



=1 a=l 



where rj^ is | rj^ - r^^ | and a—0 indicates the oxygen atom[?] . For periodic systems using 
the Ewald sum, the energy is 



'^m.olec 

E= E E|46 

i=l j<i 



a 



12 



a 



JiO,jO, 

I An 



JiOJO, 



3 3 

E E QiaQj73erfc(fi;ria,i/3)/riaj73 

a=l/3=l 



JG-Ti, 



i oc 



E E 

i=l a=l 



+ E E ( + 9 E ~ ^^^('*^i",i/3)/riQ,i/3) ) ^ (2-17) 

^ 13=1 
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where k is a screening parameter, G is a recripocal lattice vector of the periodic simulation 
cells, erf(x) is the error function, erfc(x) is the complementary error function, and L is the 
side length of the primary simulation box[?]. 

There are three independent electrostatic parameters, Xo Co, and (h, since the en- 

ergy is dependent on the difference of the atomic electronegativities and the Slater exponents 
describe J(r). We have adjusted these three parameters plus the two Lennard- Jones param- 
eters to obtain the correct gas-phase dipole moments and to optimize the energy, pressure, 
and pair correlation functions of the liquid. The parameters are given in Table |^. In order 
to implement the fluc-q procedure for the rigid bond length and rigid angle potentials used 
here, the Coulomb overlap integral needs to be evaluated only at r=0 and at the intramolec- 
ular bond lengths, tqh and thh- These values are also given on Table |l[ This optimization 
procedure does not uniquely define a set of parameters and those listed on Table |l] are one 
possible good choice which leads to improved water properties relative to the fixed charge 
models, as discussed in the next section. The electrostatic parameters, Xo ~ Xn, Jrh and 
Jqo are within the range of values used in previous EE models [?, ?, ?, ?, ?]. 

3 Numerical Method 

The molecular dynamics simulations were performed on the Connection Machine CM-5 with 
256 molecules. Periodic boundary conditions were imposed, using the Ewald sum for the 
long-ranged electrostatic potentials. The screening parameter, k, was set to 5/L and 256 
reciprocal lattice vectors were used in the Fourier space sum. A time step of 1 femtosecond 
and the SHAKE algorithm for enforcing bond constraints were used[?]. The data reported 
in the next section is from 20 separate 50 picosecond runs for each FQ model. MD is 
implemented on the parallel architecture of the CM-5 by arranging data on a 2-dimensional 
grid, so that each virtual node contains all the information (namely Tiajf3 and Qia-Qjfs) 
for the interactions between molecules i and j[?]. This is done as follows: 
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1. For a=l to a=^atoms, T^ia is placed on the first row of an N^oiec x N^oiec array (Aa(i,l)=ria) 
and tlie data is tlien spread tlirougli eacli row of tlie array using the parallel copy op- 
eration (this step takes 3% of one time step). 

2. Similarly, is placed on the first column of a different matrix (BQ,(l,i)=rjQ) and 
spread though each column(3%). 

3. Nearest image rjcj/? values can now be computed from AQ,(i,j) and B/3(i,j) without any 
communication between different virtual nodes (14%). 

4. A similar "spread-spread" algorithm is used to place Qj^ and Qj^ on each virtual node 
hi (3%) and QjQ,- is calculated (1%). 

5. All pairwise forces and potential energy terms can be calculated on each virtual node 
without communication (the Lennard- Jones term takes about 3%, the real space Ewald 
term takes 33%, and the Coulomb self-term takes 3% of a time step). The Fourier part 
of the Ewald sum (see Eq. |2.17] ) does not involve pair terms and is calculated separately 
from the N^^^g^ process (21%). 

6. The total energy and the force on atom i from all other atoms j are found from a 
parallel add operation across the virtual nodes (10%). 

7. The positions and charges are propagated using Eqs. |2]^ and|2]^(0%), bond constraints 
are enforced using SHAKE (4%), and the steps are repeated. 

On a 16 processor CM-5, we found performances of about 0.4 seconds/time step, which is 
about a factor of 10 faster than a comparable program on a IBM 580[?]. The CPU required 
for the fluctuating charge model is only a factor of 1.10 larger than for the corresponding 
fixed-charge model. 
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4 Results 



The properties of the fixed-charge and fiuc-q models for the water monomer, water dimer, 
and hquid water are hsted in Tables H and |^. Our reported error bars represent two standard 
deviation error estimates. Table |^ lists the results for the TIP4P-FQ model, comparing to 
other TIP4P geometry models: TIP4P[?], the Watanabe-Klein (WK) model [?], and the 
dipole polarizable SRWK (SRWK-P) model[?]. The SRWK model has a slightly different 
geometry than TIP4P, the length roM is 0.26A rather than 0.15A. Table ^ lists the results for 
some of the many models with an SPG geometry, including SPC[?], SPC/E [?], polarizable 
SPG (PSPC)[?], and fiexible charge SPG (SPG-FQ). The TIP4P, WK, SPG, and SPG/E 
models are all non-polarizable fixed-charge models, the WK and SPG/E models include 
a correction for the polarization energy, discused below. The WK model also includes a 
correction for the quantum librational energy. 

For the monomer, the electrostatic parameters are chosen to give the correct gas phase 
dipole moment. Another property of the isolated molecule is the dipole polarizability tensor, 
a, defined by 

^ind = • E (4.1) 

where fi'^'^'^is the dipole moment induced by the external electric field, E. ^^"•''can be de- 
termined by adding a term — ^t ■ E to Eq. |2.2| for the monomer and minimizing the total 
potential energy with respect to the charge. If the plane of the molecule is in the zy plane 
and the dipole (G2-axis) is along the z-direction, then it is found that 



2/2 



a. - 



ttxx = (4.2) 

where imh is the z-component of the distance from the negative charge site (the M site for 
TIP4P, the oxygen site for SPG) to the positive charge sites and axx is zero since all the 
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charges are in the zy plane and no dipole induction is possible out of plane. Experimentally, 
a is almost isotropic[?], so the lack of polarizability in the x-direction is clearly a deficiency 
in fluc-q models, but corrections are possible [?]. 

The dimer properties listed in Table ^ and |^ are the energy of the minimum energy 
configuration and the oxygen-oxygen distance of this configuration. Pair potentials, such 
as SPC and TIP4P, are parametrized to give the measured liquid state energies and radial 
distribution functions. It is known that this parametrization of these pairwise interaction 
models overestimates the gas phase water dimer energy. The fluctuating charge potentials 
predict an oxygen-oxygen separation closer to the experimental value but underestimate the 
dimer energy. 

We have calculated both static and dynamical properties of liquid state water at a tem- 
perature T = 298K and p = Ig/cm^. The error bars represent two standard deviations. The 
parameters for both fiuc-q models are chosen to give a binding energy of -9.9 kcal/mole. This 
energy, unlike the fixed-charge potential energies, includes the self polarization contribution 
arising from the difference in the internal energy given by Eq. for the liquid state charges 
and the gas phase charges, 

AEseif pol = E [(E,(QL^«)) - E„(Qr)" . 

a=l 

This self polarization energy is the difference between the self-energy in the liquid phase 
and the gas phase Coulombic energy. The average self polarization energy is 5.7 kcal/mole 
for TIP4P-FQ and 7.6 kcal/mole for SPC-FQ, which represents a large contribution to the 
total energy. The dipole polarizable model of Sprik and Klein has a polarization energy 
of 5.9 kcal/mole. It has been noted that because the polarization energy is ignored in the 
parametrization of the SPC and TIP4P interaction potentials, these models underestimate 
the attractive pair interactions in water [?]. One simple correction for this is to subtract 
an estimate of the polarization energy, (fiug — figas)'^ /^a, from the experimentally measured 
binding energy and to use the result to parameterize the pair potential. Here pug and figas are 
respectively the liquid state and gas phase permanent dipole moments used in the models. 
From this, it follows that the strength of the pairwise interaction must be increased to give 
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the correct energy of —10 kcal/mole. This has been done for the SPC model (giving the 
SPE/E model [?]) and TIP4P (giving the the WK model [?]), both of which have increased 
atomic charges and reduced diffusion constants compared to SPC or TIP4P. 

The pair correlation functions give detailed information about the structure of the liquid. 
The TIP4P and TIP4P-FQ oxygen-oxygen radial distribution functions, goo(r), are shown 
in Fig. and are compared to the neutron diffraction results of Soper and Phillips[?]. Recent 
experiments of Soper and Turner indicate that there is a large experimental uncertainty in the 
peak heights of the pair correlation functions, perhaps due to the use of different methods 
for removing the contribution from self or single atom scattering [?]. The peak positions 
show much less uncertainty and therefore provide more reliable points for comparison. The 
goo(r) of the flexible charge model has a first peak at a larger r than the fixed charge 
model and shows more long-ranged ordering due to the increased charges. The oxygen- 
hydrogen (Fig. ^ and hydrogen-hydrogen radial distribution function (Fig. ^ for TIP4P- 
FQ and TIP4P models are also shown. The radial distribution functions for the SPC-FQ 
potential are shown in Figs. and 0. Again, the flexible charge goo(r) shows more long- 
range correlation than the corresponding fixed charge model, but, in general, the SPC-FQ 
model does not give as accurate pair correlation functions as the TIP4P-FQ model. The 
number of nearest neighbors can be determined by integrating goo(r) over the first peak. 
Experimentally, the first minimum occurs at 3.3 A and using this value as the limit of the 
first peak gives essentially the same coordination number for all models: 4.2 (SPC-FQ), 4.3 
(TIP4P), and 4.4 (TIP4P-FQ, experiment and SPC). 

The average dipole moment, (/i), is increased for both fluctuating charge models over the 
corresponding fixed-charge models. The value of the liquid state dipole moment is not known 
experimentally. The experimental value for ice is 2.6 D[?]. Theoretical studies which use 
potentials which have the correct dipole polarizability and quadrupole moments find a dipole 
moment of 2.5 D [?] and 2.45 to 2.7 D, depending on the details of the non-electrostatic part 
of the potential [?]. Additionally, it has been observed that the dependence of the dielectric 
constant on the dipole moment is such that to have a dielectric constant close to 80, the 
potential must have a dipole moment in the range of 2.3 to 2.6 D[?]. The distributions of 
the dipole moment are shown in Figure The full width at half maximum is 0.42 for the 
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SPC-FQ and 0.49 for TIP4P-FQ. The SPC geometry, with greater distances between the 
charge sites, is more polarizable and therefore the distribution of the \fi\ is broader. In the 
fluctuating charge models the instantaneous dipole moment does not always lie along the C2 
direction (what we have previously defined as the z-axis). The average of the component of 
the dipole moment along the C2 axis, fi^, for the TIP4P-FQ model is 2.59 D, which when 
compared to the total dipole moment of 2.62 D indicates that there are small fluctuations 
of the dipole moment away from the C2 axis. For the SPC-FQ model, fi^ is 2.81 D and the 
total dipole moment is 2.83 D. 

The static dielectric constant, eo, for the FQ potentials, calculated from the fluctuations 
in the total dipole of the central simulation box, M, according to[?] 



is 79±8 for TIP4P-FQ and 116±18 for SPC-FQ. Eq. |4.3| was evaluated from 1 nanosecond 
run. The dielectric constant of the TIP4P-FQ model is in good agreement with experiment, 
which is consistent with earlier findings that models with a dipole moment of 2.6 D have 
a dielectric constant near 80[?, ?, ?]. The SPC-FQ model is in poorer agreement due to 
the fact that the liquid state dipole moment is large. A large dipole moment[?, ?, ?] and 
dielectric constant [?] is also seen in dipole polarizable SPC models. From this perspective 
the TIP4P geometry is better than SPC. 

The fluctuating charge models have a dielectric response characterized by an infinite 
frequency dielectric constant, e^o- Neumann and Steinhauser have derived expressions for 
calculating eoo for dipole polarizable models[?]. Expressions for eoo for fluctuating charge 
models are similar. The Coulomb energy can be written as 



where Q is a 2 Nmoiec dimensional vector containing the hydrogen charges (the oxygen charges 
are eliminated from the equation using the charge neutrality constraint), I is the identity 




(4.3) 



u(Q) = (x?i-x;^)Q-i + ^Q-J-Q 



(4.4) 
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vector, and J is a (2Nmoiec x 2Nmoiec) matrix. The elements of J are given by 

With the Ewald sum — which we are using in the present calculations — for the long range 



Coulomb interaction, each 1/r and Jq^ interaction in Eq. ^]5| is replaced as follows, [?] 

L (or 

Ja,l3(jia,il3) =^ J a, (3 10,1/3) — ^'^^{t^'^ia,ip) / '^ia,ip + TT X! ^^'^ COs(G ■ Fja^j^). 

The minimum energy charges are given by 

Q = -(x?i-X°o)J-^-I (4.6) 

where is the inverse of J. In the presence of a spatially homogeneous external electric 
field, E, the charges are 

Q = qo + J"^ ■ 5r ■ E (4.7) 
where Q° are the charges in the absence of the field and (5rja=ria — Yio- The energy is 

U = U°-M-E--E-A-E (4.8) 

where U° is the energy in the absence of the field and A is the polarizability of the system, 
given by 

Kajp = Sri^ ■ SrjpJ-^jp. (4.9) 

Following Ref. [?], e^o is 

A/jr -^molcc -^molcc 3 3 

e°o = i + w E E EE(A.../3). (4.10) 

i=l j=l a=2/3=2 

Eq. [4.10| was evaluated every 50 picoseconds in the course of the simulations, which is a 
frequent enough sampling rate to provide a precise estimate of {Kiajp) . The value obtained 
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for eoo from both fluctuating charge models is about 1.6, close to the experimental value 
of 1.79 [?]. 600 is underestimated because the perpendicular polarizability, Oxx, is zero. To 
leading order in a, and thus eoo-1, the total system polarizability, A, is proportional to 
Tr(Q:)[?]. It then follows that to order a, 

6oo(FQ)-l ^ Tr[a(FQ)] 



€00 (exact) — 1 Tr [a (exact)] 

which is true for the data from Tables ^ and |^. Therefore, = produces an eoo smaller 
than the experimental value. 

The frequency dependent dielectric constant can be calculated from[?] 

^^^ = 1-^.; A. [0(^)1 (4.12) 

where denotes the Laplace operator and (f){t) is the normalized time autocorrelation 
function of the system's total dipole (M=X] t-ti), 

0(t) = (M(t)-M(O))/(M2). (4.13) 

0(t) has a short-time oscillatory part, due to librational motions of the hydrogen atoms. At 
long times, (f){t) decays exponentially and the decay constant is the Debye relaxation time, 
td (see Fig. In order to perform the Laplace transform to get €{uj), we set 0(t)=Aexp[- 
t/rn] for times longer than 0.5 picoseconds. The parameters A and tq were chosen to give a 
smooth interpolation between the calculated 0(t) and the exponential fit (see Fig. P). The 



frequency dependent dielectric constant for the TIP4P-FQ model is shown in Fig. |T0[ The 
agreement with the experimental results [?, ?, ?] is very good. The close agreement in the 
low-frequency microwave range is due to the fact that the model gives accurate values of eo 
and td. The features at frequencies higher than 300 psec~^ are due to bond stretches and 
bends and so are not present in the rigid geometry models used here. The highest frequency 
feature given by the TIP4P-FQ model is the librational mode, which shows a peak in e" 
at 130 psec~^. The experimental peak is at 90 psec~^[?]. Another notable feature is at 25 
psec~^ which has been interpreted as a translational vibration of a water molecule in its cage 
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of nearest neighbors[?]. This feature is not present in the spectrum for non-polarizable water 
models such as TIP4P[?] or MCY[?]. As argued by Neumann, this translational motion will 
not change the system's dipole moment much for non-polarizable models, but for polarizable 
models, the translation motion will induce a change in the dipole moment [?]. Therefore this 
feature can only be seen in polarizable models. The fluctuating charge models support that 
argument. The frequency dependent dielectric constant for SPC-FQ is shown on Fig. [TTi 
The agreement with experiment is not as good as the TIP4P-FQ model, primarily because 
the static dielectric constant is overestimated. The librational peak is at 160 psec~^ for the 
SPC-FQ model. 

Lastly, we examine the dynamical properties of the fixed and fiuc-q model potentials. 
In general, the flexible charge models have slower translational and rotational time-scales 
than the fixed-charge models, primarily due to the stronger electrostatic interactions from 
the higher charges. The translational diffusion constant, D, is determined from the Einstein 
relation 

D= liml(|rf^(t)-rf^(0)p) (4.14) 

where rp^(t) is the position of the center-of-mass of molecule i at time t. The diffusion 
constants for the flexible charge models are smaller than the fixed-charge models and closer 
to the experimental value (see Tables ^ and ^). Rotational time constants are calculated 
from 

cut) = {Piietit) ■ etm) (4.15) 

where P; is a Legendre polynomial and e" are unit vectors along molecule ?'s principle axis 
of rotation[?]. Rotations around the axis connecting the hydrogen atoms (the y-axis) can 
be measured by proton NMR. The zero frequency component of the Fourier transform of 
C2(t) gives the NMR relaxation time, t^va/r- The correlation function C^t) has a long range 
exponential decay, given by exp(-t/r2), and a short range, oscillatory part (see Fig. [T^) . 
The zero frequency part of the Fourier transform of C^t) is, to a good approximation, given 
by A2 r|, since the short range part will not contribute much to the transform. The flexible 
charge models exhibit slower reorientational dynamics than the fixed-charge models and are 
in closer agreement with the experimental value [?]. 
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5 Conclusions 



The fluctuating charge (fluc-q) water models, using either a TIP4P or SPC geometry, were 
shown to give important improvements over fixed charge models. The electronic properties 
of the fiuc-q models are such that in the gas phase they give the correct dipole moment (this 
is by construction) and in the liquid phase, the dielectric properties are well reproduced 
for a range of frequencies (see Figs. |l^ and |Tl|). The dielectric properties are best for the 
TIP4P-FQ model which gives a static dielectric constant, eo, an infinite frequency dielectric 
constant, eoo, and a Debye relaxation time, r/j, close to the experimental values. The SPC- 
FQ model overestimates eo, but and td are accurate. The fiuc-q models also show (at 
25 psec~^) a feature in the dielectric spectrum originating from translational motion of a 
water molecule in the cage of its neighbors. This is a feature which fixed charge models do 
not show[?, ?]. Translations will have a large effect on the system's dipole moment only 
for polarizable models, so this feature of e{ui) is an indication of the coupling between the 
electronic and nuclear degrees-of-freedom. 

In addition, the fiuc-q models give good estimates for the liquid-state radial distribution 
functions (see Figs. 0- 1^ and dynamical properties such as the diffusion constant and tnmr 
(see Tables |^ and Since the charges are not fixed to values which represent mean field 
values for a particular single-component phase, this method should be transferable to studies 
of heterogeneous systems in which deviations from the mean field charge values should be 
greater than in pure systems[?]. The fiuc-q method assigns two parameters to each element 
corresponding to the two terms in the energy expansion (Eq. |2.1j ). The first order term 
is the Mullikan electronegativity, and the second order term is determined by a Slater 
exponent, (, through the Coulomb overlap integral (Eq. |2.11|) . Extensions to more complex 
molecules would require additional terms for other elements, perhaps to be taken from other 
electronegativity equalization schemes [?, ?, ?, ?, ?]. 

All in all, the fiuc-q water models are as successful as the best of the dipole polarizable 
models, such as SRWK-P[?], RER[?], and RPOL[?]. The dipole polarizable models intro- 
duce an interaction (the dipole-dipole interaction) not present in fixed-charge models, which, 
together with solving for the induced dipole moments, increases the CPU time by about a 
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factor of 2 over fixed charge models[?, ?]. The fluc-q method introduces no new interactions 
and the propagation of the charges using extended Lagrangian methods increases the com- 
putational cost by only a modest amount (about 1.1). The dynamical fluctuating charge 
method is therefore a tractable means for the inclusion of charge polarization effects and will 
be useful for the study of large systems, such as proteins in aqueous environments. 
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Table 1: Potential parameters for the fixed-charge potentials, SPC and TIP4P, and the 
flexible charge SPC and TIP4P (SPC-FQ, TIP4P-FQ) models. The last four terms, Jab, are 
determined by and Co and so are not independent parameters. 





SPC" 


TIP IP'' 


SPC-FQ 


TIPIP-FQ 


e (kcal/mole) 


0.1554 


0.1550 


0.2941 


0.2862 


-(A) 


3.166 


3.154 


3.176 


3.159 


dnoH (degrees) 


109.47 


104.52 


109.47 


104.52 


ro// (A) 


1.0 


0.9572 


1.0 


0.9572 


roM (A) 


0.0 


0.15 


0.0 


0.15 


Qi/ (e) 


0.41 


0.52 






Xo - Xh (kcal/mole e) 






73.33 


68.49 


Co (ao^) 






1.61 


1.63 


Ch (ao ^) 






1.00 


0.90 


J^o (kcal/(mole e^)) 






367.0 


371.6 


J%fj (kcal/(mole e^)) 






392.2 


353.0 


Jof/(roH) (kcal/(mole e^)) 






276.0 


286.4 


Jhh{^hh) (kcal/(mole e^); 


) 




196.0 


203.6 



a) Ref. [?], b) Ref. [?]. 
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Table 2: Properties for potentials with the TIP4P geometry: the fixed-charge TIP4P and 
Watanabe-Klein (WK) models, the SRWK dipole polarizable model (SRWK-P), and the 
flexible charge model (TIP4P-FQ). Properties listed are the the gas-phase dipole moment, 
the dipole polarizabilities, an (the y and z directions lie in the plane of the molecule, with 
the z-axis along the C2 axis), the energy of the dimer in its minimum energy configuration, 
the distance between oxygen atoms for the minimum dimer configuration, and properties of 
the liquid as indicated. 





TIP4P" 


WK'' 


SRWK-P'^ 


TIP4P-FQ 


experimental 


Gas-phase dipole moment (Debye) 


2.18 


2.60 


1.85 


1.85 


1.85^* 


(A3) 








1.44 


0.82 


1.468±0.003^ 


ayy (A3) 








1.44 


2.55 


1.528±0.013^ 


(A3) 








1.44 





1.415±0.013^ 


Dimer energy (kcal/mole) 


-6.3 






-4.5 


-5.4±0.7^ 


Dimer 0-0 length (A) 


2.75 






2.92 


2.98-'^ 


Liquid state properties (T=298 K, p—l.Q g/cm3) 








Energy (kcal/mole) 


-10.1" 


-10.2^ 


-11. P 


-9.9 


-9.9" 


Pressure (kbar) 


0.0^ 


O.l'' 


0.6^ 


-0.16±0.03 


0.0 


Dipole moment (Debye) 


2.18 


2.60'' 


2.63" 


2.62 




eo 


53±2» 


80±8'' 


86±10" 


79±8 


78'^ 




1 


1 




1.592±0.003 


1.79'^ 


Diffusion constant (10~^ m^/s) 


3.6±0.2^ 


1.1±0.3'' 


2.4±0.3'^ 


1.9±0.1 


2.30^ 


tnmr (ps) 


1.4±0.2* 


3.8±0.3'' 




2.1±0.1 


2.V 


td (ps) 


7±2'' 


22±4'' 




8±2 


8.27±0.02^ 


l) Ref. [?], b) Ref. [?], c) Rcf. [?], 


d) Ref. [?], 


e) Ref. [?; 


1, f) Ref. [?: 


1, g) Ref. [?], h) Ref. 



[?],i)Ref. [?],j)Ref. [?], k) Ref. [?] 
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Table 3: Properties for potentials with the SPC geometry: the fixed-charge SPC and SPC/E 
models, dipole polarizable model (PSPC), and the flexible charge model (SPC-FQ). Proper- 
ties listed are the the gas-phase dipole moment, the dipole polarizabilities, an (the y and z 
directions lie in the plane of the molecule, with the z-axis along the C2 axis), the energy of 
the dimer in its minimum energy conflguration, the distance between oxygen atoms for the 
minimum dimer conflguration, and properties of the liquid as indicated. 





SPC* 


SPC/E'' 


PSPC^ 


SPC-FQ 


experimental 


Gas-phase dipole moment (Debye) 


2.27 


2.35 


1.85 


1.85 


1.85^* 


a.. (A3) 








1.44 


1.02 


1.468±0.003'= 


(A3) 








1.44 


2.26 


1.528±0.013^ 


axx (A^) 








1.44 





1.415±0.013^ 


Dimer energy (kcal/mole) 


-6.7 






-4.4 


-5.4±0.7-'^ 


Dimer 0-0 length (A) 


2.75 






2.94 


2.98^ 


Liquid state properties (T=298 K, p 


=1.0 g/cm 










Energy (kcal/mole) 


-10.0^? 


-Q.Q'* 


-9.F 


-9.9 


-9.9' 


Pressure (kbar) 


0.3» 


-0.08i0.04'* 




0.03±0.05 


0.0 


Dipole moment (Debye) 


2.27 


2.35 


2.9" 


2.83 






68±7^' 


67±10'* 




116±18 


78^^ 




1 


1 




1.606±0.002 


1.79^' 


Difi'usion constant (10^^ m^/s) 


3.3±0.2^' 


2.4±0.4'* 


2.0±0.2'^ 


1.7±0.1 


2.30' 


tnmr (ps) 


l.l±0.2s 


l.QiO.l'* 




2.2±0.1 


2.1"* 


td (ps) 


ll±2s 


10±3'* 




9±3 


8.27±0.02" 



a) Ref. [?], b) Ref. [?], c) Ref. [?], d) Ref. [?], e) Ref. [?], f) Ref. [?], g) Ref. [?], h) Ref. 
[?], i) Ref. [?], j) Ref. [?], k) Ref. [?], 1) Ref. [?], m) Ref. [?], n) Ref. [?] 
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Figure 1: Coulomb interaction for H-H pairs, comparing the intramolecular Coulomb overlap 
interaction (solid line) with the intermolecular pure Coulomb, 1/r, interaction (dotted line), 
in kcal/mole/e^. 



Figure 2: Oxygen-oxygen radial distribution function for the TIP4P-FQ (solid line) and 
TIP4P (dotted line) potentials, comparing to the neutron diffraction results of Soper and 
Phillips (dashed hne). 



Figure 3: Oxygen-hydrogen radial distribution function for the TIP4P-FQ (solid line) and 
TIP4P (dotted line) potentials, comparing to the neutron diffraction results of Soper and 
Phillips (dashed line). 
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Figure 4: Hydrogen- hydrogen radial distribution function for the TIP4P-FQ (sohd Une) and 
TIP4P (dotted Une) potentials, comparing to the neutron diffraction results of Soper and 
Phillips (dashed line). 



Figure 5: Oxygen-oxygen radial distribution function for the SPC-FQ (sohd hue) and SPC 
(dotted line) potentials, comparing to the neutron diffraction results of Soper and Phillips 
(dashed line). 



Figure 6: Oxygen- hydrogen radial distribution function for the SPC-FQ (sohd hne) and SPC 
(dotted line) potentials, comparing to the neutron diffraction results of Soper and Phillips 
(dashed hne). 



Figure 7: Hydrogen-hydrogen radial distribution function for the SPC-FQ (solid line) and 
SPC (dotted line) potentials, comparing to the neutron diffraction results of Soper and 
Phillips (dashed hne). 



Figure 8: Distribution of the absolute value of the dipole moments of the TIP4P-FQ (solid 
line) and SPC-FQ (dashed line) liquids. 
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Figure 9: Time autocorrelation function for the system dipole for the SPC-FQ (top Une) and 
TIP4P-FQ (bottom Une) models and the exponential long-time approximation (dash lines), 
shown on semilog axes. 



Figure 10: Real (top) and imaginary (bottom) parts of the frequency dependent dielectric 
constant for the TIP4P-FQ model (sohd hues), compared to experiment (dotted hues). 



Figure 11: Real (top) and imaginary (bottom) parts of the frequency dependent dielectric 
constant for the SPC-FQ model (sohd lines) , compared to experiment (dotted hues) . 



Figure 12: Rotational correlation function for TIP4P-FQ (sohd line) and SPC-FQ (dotted 
line) . 
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